library(SDMTools)#load necessary libraries
work.dir = 'C:/Users/alice hughes/Documents/Pakistan_models/2050A1B/overlap/'; setwd(work.dir) #define & set the working directory

###01.create a vector of species names
files=list.files(pattern='_overlap.asc') #get a list of files in your working directory that include 'asc'

richness=NULL
for (tfile in files) {cat(tfile, '\n')
	overlap=read.asc.gz(tfile)
	overlap[which(overlap<3)]=0  #set values less than 3 to zero
	overlap[which(overlap==3)]=1  #set values equal to 3 (where current and future overlap) to 1
	
	if (is.null(richness)) {richness = overlap;  } #define the output & set everything = 0
	richness = richness + overlap
	}
dev.off()

###create binary ascii and store in memory

write.asc.gz(richness, "overlap_bat_richness.asc") #write out the richness data
base.asc=richness
base.asc[which(is.finite(base.asc))]=0
richness[which(richness==0)]=NA
max.richness = max(richness, na.rm=TRUE)
min.richness=min(richness,na.rm=TRUE)
#threshold=1
#bins = seq(0,1,length=101); bins = cut(threshold,bins,labels=FALSE) # get the threshold bin for cols
#cols = c(rep('gray86',bins),colorRampPalette(c("wheat","tan","forestgreen","darkgreen"))(100)[bins:100])

cols = colorRampPalette(c("tan","forestgreen"))(100)

 #get the max richness value
png('overlap_bat_richness.png', width=nrow(richness),height=ncol(richness), units='px', res=300, pointsize=5, bg='white') #start the plot
	par(mar=c(0,0,0,0)+0.1) #remove any plot margins
	image(base.asc,ann=FALSE,axes=FALSE, col='gray86')
	image(richness, ann=FALSE,axes=FALSE,col=cols, zlim=c(min.richness,max.richness),add=TRUE) #plot the richness
	legend.gradient(legend.pnts,cols=cols,limits=c(0,max.richness), title='Richness', cex=2)
dev.off()
